Mon. Not. R. Astron. Soc. 000, 1-11 (2012) Printed 11 October 2012 (MN lAT^X style file v2.2) 



Resonances and bifurcations in systems with elliptical 
equipotentials 

Antonella Marchesiello 1 * and Giuseppe Pucacco 2 ' 3 f 

1 Dipartimento di Scienze di Base e Applicate per I'Ingegneria - Universitd di Roma "la Sapienza" 

2 Physics Department, University of Rome "Tor Vergata", 100133 Rome 
3 INFN, Sezione Roma Tor Vergata 



Received 19/09/2012. Accepted 08/10/2012. 



ABSTRACT 

We present a general analysis of the orbit structure of 2D potentials with self- 
similar elliptical equipotentials by applying the method of Lie transform normaliza- 
tion. We study the most relevant resonances and related bifurcations. We find that 
the 1:1 resonance is associated only to the appearance of the loops and leads to the 
destabilization of either one or the other normal modes, depending on the ellipticity 
of equipotentials. Inclined orbits are never present and may appear only when the 
equipotentials arc heavily deformed. The 1:2 resonance determines the appearance of 
bananas and anti-banana orbits: the first family is stable and always appears at a 
lower energy than the second, which is unstable. The bifurcation sequence also pro- 
duces the variations in the stability character of the major axis orbit and is modified 
only by very large deformations of the equipotentials. Higher-order resonances appear 
at intermediate or higher energies and can be described with good accuracy. 
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1 INTRODUCTION 

In investigating the orbit structure of a galactic potential, we 
are often interested in some particular feature of its general 
layout. We may mention: the birth and/or disappearance of 
specific orbit families; their stability nature; the phase-space 
fraction occupied by invariant tori around stable periodic or- 
bits, etc. (Binney & Tremaine 2008). In integrable systems 
these features are uniquely determined by the integrals of 
motion (de Zeeuw 1985b): only a limited number of orbit 
families exist and their possible bifurcations occur at iso- 
lated critical values of the conserved functions giving the 
integrals. 

On the other hand, the dynamics of generic systems are 
not integrable. There are several bifurcations with a prolif- 
eration of periodic orbit families and sooner or later a tran- 
sition to a stochastic behavior. Stochasticity, if not limited 
to small regions of phase space, leads to chaos (Contopoulos 
2004) . However non-integrable dynamics do not prevent reg- 
ular behavior: significant parts of phase space can be layered 
with invariant surfaces and in many instances a generic sys- 
tem as a whole can be quite similar to an integrable system 
(Henon & Heiles 1964). In these circumstances perturbation 
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approaches can be devised to describe the features of the 
system (Gustavson 1966). 

A powerful perturbation method is that based on 
Hamiltonian Normal Forms (Boccaletti & Pucacco 1999). 
Typically, the application of this method is based on three 
steps: 

1. Construction of a new Hamiltonian (the 'normal form') 
by means of a canonical transformation suitable to capture 
a peculiarity of the system under study. 

2. Use of the normal form to investigate in the simplest way 
the particular feature of the system we are interested in. 

3. Inversion of the transformation to describe this feature 
in terms of the original parameter of the system allowing, if 
possible, the comparison with observational data. 

In some cases, steps 2. and 3. can be reversed but, usu- 
ally, working with the normal form in the normalization co- 
ordinates is easier and/or more effective. 

We recall that coping with non-integrable dynamics 
through perturbation theory often means to try to com- 
pute non-existing quantities. To clarify this seemingly ab- 
surd statement, we can say that in any perturbative ap- 
proach dynamical quantities are expressed as series in some 
(small) parameter. Physical properties of the original system 
(e.g. the gravitational potential) give convergent series in 
a suitable neighborhood. The normalization procedure pro- 
vides expansion series for quantities approximating phase- 
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space invariants. As a matter of fact, these series do not con- 
verge: if they would, this should imply that global invariants 
existed and that the dynamics should be everywhere regu- 
lar, which is not. However, due to the asymptotic nature of 
these series (Sanders et al. 2007), their truncations can be 
used to approximate local invariant surfaces responsible for 
regularity. The change in their geometric properties caused 
by the variation of physical parameters allows us to locate 
bifurcations of new orbit families. A remarkable fact is that 
often the neighborhood of asymptotic semi-convergence is 
large enough to provide reliable information in a wide range 
of physical parameters and several features of non-integrable 
dynamics can be accurately predicted (Pucacco et al. 2008) . 
Technical issues do not always lead to a straightforward ap- 
plication of the method: just to mention a few, we recall 
that among the many features characterizing the real sys- 
tem, the normal form is able to describe only a limited subset 
of them, typically in the neighborhood of a given resonance 
(Bclmonte et al. 2007). Another issue is that there is no 
definite strategy to predict the best truncation of the series 
expansions; moreover, it can be not easy to re-express them 
in terms of observables. Generic cases (colloquially, models 
with 'many' parameters) are cumbersome: mathematicians 
have therefore introduced simplifying techniques (singular- 
ity theory, catastrophe theory, etc.). However, in spite of 
their power and elegance, they are even more difficult to 
use so that, in applications, they are still not so useful. We 
prefer to stay on 'standard' methods. 

Aim of the present paper is to offer a well defined set- 
ting in which many of the technical issues listed above are 
addressed and solved. We will see how to construct normal 
forms for the dominating resonances, how to use higher- 
order expansions to predict bifurcation thresholds and sta- 
bility transitions and will show circumstances in which these 
quantities can be computed in a 'large' range of parameters. 
Resonance between two non-linear oscillations is the source 
of non-trivial dynamics (Contopoulos 1963; Contopoulos & 
Moutsoulas 1966; Verhulst 1979; Binney 1981). de Zeeuw 
& Merritt (1983) made a general analysis of the symmetric 
1:1 resonance with the averaging method. The method of 
normal form is more flexible in treating generic resonances 
requiring higher-order computations and, when applied to 
the same models, the results are identical to those of the 
averaging method (Marchesiello & Pucacco 2011). However, 
when a comparison with numerical results (see e.g. Miralda- 
Escude & Schwarzschild (1989)) requires precise predictions, 
higher-order computations are necessary (Belmonte et al. 
2007). The setting in which we work is that of potentials 
with similar concentric ellipsoidal equipotentials. To shed 
light on the methods and to limit the algebraic complica- 
tions we limit the treatment to 2-DOF non-rotating systems. 
Explicit formulas for the bifurcation thresholds of the main 
periodic orbits are computed in terms of the energy for a 
family of models with two shape parameters. Additional el- 
lipsoidal symmetry-breaking perturbations are included. We 
also discuss the relation with other issues like Stackel fits to 
separable systems, surfaces of section, order and chaos, etc. 
In particular, we stress how the asymptotic nature of the 
method allows us to make reliable predictions in a domain 
much larger than expected on the basis of standard pertur- 
bation arguments. 

The plane of the paper is as follows: in Section 2 we in- 



troduce the procedure to construct the approximating inte- 
grable systems by recalling the method of the Lie transform 
(Gerhard & Saha 1991; Giorgilli 2002); in Sections 3 and 
4 we apply this approach to investigate the main aspects 
of the dynamics in a symmetry plane of a triaxial ellipsoid 
(Belmonte et al. 2007) , obtaining first-order estimates of the 
bifurcation thresholds of the 1:1 and 1:2 periodic orbits; in 
section 5 we analyze higher-order cases; in Section 6 we dis- 
cuss further developments and hints for other applications 
and in Section 7 we conclude. 



2 RESONANT HAMILTONIAN NORMAL 
FORMS 

Let us consider a two degree of freedom system with a 
smooth potential with an absolute minimum in the origin, 
symmetric under reflection with respect to both coordinate 
axes. The Hamiltonian is given by 



H{w) = \{ P l+p 2 y ) + V {N) {x\y 2 ) 



(1) 



where with w we collectively denote the phase-space vari- 
ables and we assume that the potential can be expanded as 
a truncated power series 



vW(*V) = $>*(;nV) 



where 



fc+i 



v,( 1 2 , J 2 ) = Ecw-^ 2, / (M 



(2) 



(•3) 



3=0 



mined by the problem under study. 

In particular, we are interested in a fairly general class 
of potentials with self-similar elliptical equipotentials of the 
form 



V(x,y;q,a) 



i ( 1 + x 2 + \ 



a /a 



log ( 1 + X 



2 , y- 



< a < 2 
a = 0. 



(4) 



The ellipticity of the equipotentials is determined by the 
parameter q: for short, we will speak of an 'oblate' figure 
when q < 1 and a 'prolate' figure when q > 1. The profile 
parameter a determines the behavior at large radius. 

The family of potentials (4) can be expanded in a series 
of the form (3), or more simply of the form 



where we have introduced the 'elliptical radius' 



8(q) = Jx* + * s . 



(5) 



(6) 



With unit 'core radius' we can put Bo = 1/2 and, for the 
class (4) , the first two coefficients of the higher-order terms 
are 



2-a (2-a)(4-q) 
Bl = 8~' B2 = 48 ' 
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Another interesting case is that of the 'flattened isochrone' 
(Evans, de Zeeuw & Lynden-Bell 1990): 



B 2 = — . 

32 



(8) 



Each term in the series is given by an even power of the basic 
elliptical radius and the Hamiltonian (1) can be treated in 
a perturbative way as a non-linear oscillator system. 

To find the normal form we have first of all to 'prepare' 
the Hamiltonian. We start by introducing a small parameter 
e > and, by performing a 'blowing-up' of the phase-space 
by means of the transformation 



(0) 



we rescale the Hamiltonian (1) according to 



1 N 
H(w) = e- 2 H(w) = -(pl+ pi) + * V» (x 2 ,y 2 ). (10) 

fc=0 

With this trick we assign an order to the terms in each 
series without making an explicit reference to the extent of 
the neighborhood of the equilibrium. After a further scaling 

(11) 
(12) 



p x = ^/cJTpi, x = xi/>Ju)i, 

Py = V^2P2, y = X 2 /y/i02, 

where 



wi = V2Eh = 1, ui2 = V2Eh/q = 1/q, (13) 
the original Hamiltonian system (1) is put into the form 



H(w) = ^ S 2n H 2n , 



(14) 



where we still use w to denote phase-space variables. We 
then have 



H = - (wi(p? + x\) + to 2 (p 2 2 + x 2 2 )) 



(15) 



and H 2 j(w),j > 0, are essentially the higher (than the 
second) order terms of the potential. We are interested in 
the behavior of the system 'around' m/n resonances with 
m, n £ N: in general our frequency ratio q is an irrational 
number and the unperturbed system is non-resonant. How- 
ever the nonlinear higher-order terms produce a passage 
through resonance with interesting dynamics. To describe 
this phenomenon we introduce a 'detuning' parameter 5 
(Verhulst 1979; de Zeeuw & Merritt 1983) such that the 
frequency ratio is written as 



(16) 



uji m 

— = q = ho. 

Ul2 n 

The detuning parameter is treated as a term of order two in 
e (8 = oe ) and considered 'small'. After a further rescaling 

H = —H = 2qH (17) 

0J 2 

and noting that, in view of (16), we have 

- = 5-de H ri£ +..., (18) 

q m m l m s 

by collecting terms up to order 2N in e, we finally put the 



Hamiltonian into the form 

JV 



(19) 



where the unperturbed term (in exact resonance) is given by 

(20) 



Ho = \m{p\ + x\) + ^n(p 2 2 + x\). 

The system is now ready for a standard resonant normaliza- 
tion: it undergoes a canonical transformation to new vari- 
ables W(w), such that the new Hamiltonian is 



K(W) = z 2 " K ^ = e CG H(w), 

n=0 

where the linear differential operator 
p Cg - —r k 



(21) 



(22) 



associated to the generating function G(w), is defined by its 
action on a generic function F(w) by the Poisson bracket: 



£ G F = {F,G}. 



(23) 



To construct K starting from H is a recursive procedure 
exploiting an algorithm based on the Lie transform (Boc- 
caletti & Pucacco 1999; Giorgilli 2002). A short account 
useful for the present purpose is given in Marchesiello & 
Pucacco (2011). To proceed we have to make some decision 
about the structure the new Hamiltonian must have, that 
is we have to chose a normal form for it. We construct the 
new Hamiltonian in such a way that it admit a new integral 
of motion, that is we consider a certain function, say I(w), 
and impose that 



{K,I} = 0. 



(24) 



The usual choice (but not necessarily the only one) is that 
of assuming 

I = H = K (25) 

so that the function (20) plays the double role of fixing the 
specific form of the transformation and assuming the status 
of second integral of motion. 

Formally, a more direct way of applying this method 
is by using smarter coordinates which greatly simplify the 
procedure. A first choice is that given by the complex coor- 
dinates 



Zl 



Pi + ixi, 



22 



P2 +ix 2 , 



(26) 



leading to a normal form K(z%, z 2 ,zi, z 2 ) so that, for exam 
pie, 

1 



Ho = Ko = - {mz\z\ + nz 2 z 2 ). 



(27) 



A second useful choice is the action-angle-/iA;e variables 
J a , 6 a defined through the transformations 



Za — 2JqG 

In this way we have 
Ho = A'o = mJi + nj 2 , 
so that 



1,2. 



C 



d 



H 



+ n- 



d 



(28) 
(29) 
(30) 
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With the choice (25), condition (24) translates in the neces- 
sary condition 



C Hn K = 



(31) 



the new Hamiltonian must satisfy. Since a generic polyno- 
mial series turns out to be a Fourier series in the angles with 
coefficients depending on the actions, the typical structure 
of the resonant normal form (truncated when the first reso- 
nant term appears) is (Sanders et al. 2007) 



K = 



i,Ji+n(J 2 +e 2 8J 1 ) + ]T s 2k r {k+1) (Ji,J 2 ) + 



e 2( m+ n-i) AmnJ nj^ cos[2(nei _ m6 , 2)]j ( 32) 

where V^' are homogeneous polynomials of degree j whose 
coefficients may depend on 5 and the constant A mn {q,a) is 
the coefficient of the resonant term. It is easy to check that 
this is the most general form of a phase-space function of 
degree m + n in the actions which stays in the kernel of Ch (i 
as given by (30). In these variables, the second integral is 
given by (29) and the angles appear only in the resonant 
combination n6\ — m6 2 : for a given resonance, these two 
statements remain true for arbitrary 

N>N r = m + n-l, (33) 

where N r is the order of the resonance. New variables 
'adapted to the resonance' (Sanders et al. 2007) are intro- 
duced by means of the quasi-canonical transformation, 



(34) 



TZ= <£■ (nJi — mj 2 ), 
tp = n(n8\ — m9 2 ), 
X = /i(m0i + 716*2). 

The transformation is canonical when A = 1/fj,, but other 
choices can be convenient to simplify formulas: we will usu- 
ally choose /i — 2. Under transformation to these new vari- 
ables, the Hamiltonian can be expressed in the reduced form 

K(R,,ii;S) = uK{J 1 (TZ,8),J 2 {TZ,£),2(ne 1 -m6 2 )), (35) 

with v a scaling factor chosen to get the simplest expression 
from the quasi-canonical transformation. We obtain a family 
of 1-dof systems in the phase-plane 1Z, ip, with equations of 
motion 



■R 



dIC 
9i> • 



(36) 
(37) 

parametrized by £ that is conserved because is proportional 
to the value of the integral of motion (29). 

The dynamics of the 1-dof Hamiltonian K,(JZ,^) are 
integrable. Unfortunately, this does not necessarily implies 
that the solution of the equations of motion can be writ- 
ten explicitly. However, a quite general description of the 
phase-space structure of the system is possible if we know 
the nature of the fixed points, since these turn out to be 
the main periodic orbits of the unreduced system. In fact, 
centers (namely maxima and minima of K.) are associated 
with stable periodic orbits which parents quasi-periodic or- 
bits with essentially the same properties, whereas saddles 
of K, are associated with unstable periodic orbits. For the 
main periodic orbits of the Hamiltonian (32) , J a ,9 a are true 
action-angle variables and so the solutions to which they cor- 
respond are known. There are two types of periodic orbits 



that can be easily identified by means of the fixed points of 
the system (36,37): 

(i) The normal modes, for which one of the Jt vanishes: 
the solutions 1Z = ±8, J 2t i = are respectively the periodic 
orbit along the rc-axis and the y-axis. 

(ii) The periodic orbits in general position, namely those 
solutions characterized by fixed relations between the two 
angles, tpo = 2(n0i — m6 2 ). These are solutions of 1Z = 
(when 1Z 7^ ±£) and determine the corresponding solutions 
of 



rl>-. 



QIC 



0. 



(38) 



For all cases treated below they fall in two classes: ip = 
(to which we refer as the in phase oscillations) and ipo = ±w 
(the anti-phase oscillations). 

As a rule, normal modes exist on every surface K = 
(n/oj 2 )E, where E is the true energy. Periodic orbits in gen- 
eral position exist instead only beyond a certain threshold 
and we speak of a bifurcation ensuing from a detuned res- 
onance. The bifurcation is usually described by a series ex- 
pansion of the form 



(39) 



where the cu are coefficients depending on the resonance ra- 
tio and the parameters of the system. Eq.(39) implies that 
at exact resonance (vanishing detuning) the bifurcation is 
intrinsic in the system and that, when we move away from 
the exact ratio, the critical value E c of the threshold en- 
ergy gradually increases. We will see that already a linear 
relation given by the first order truncation provides a re- 
liable estimate of the threshold values. Actually, by using 
Hamiltonian (32), the thresholds naturally appear in terms 
of the 'distinguished' variable £: to arrive at expressions of 
the form (39) we need to disentangle the relation between £ 
and the true energy (Belmonte et al. 2007). 

By plotting curves (39) in the (q, _E)-plane we get in- 
formation for a given value of the other morphological pa- 
rameters (a in our reference cases). Each resonance corre- 
sponds to a family of periodic orbits to which it is custom- 
ary to assign the nicknames introduced by Miralda-Escude 
& Schwarzschild (1989). The nature of the critical points 
of the system (36,37) determines the stability/instability 
property of the orbit. With obvious limitations due to a 
perturbative approach we may deduce the main aspects of 
the phase-space structure. We recall that the most obvious 
limitation of the method is determined by the values of dy- 
namical and/or morphological parameters beyond which the 
dynamics are mostly chaotic. We can increase the precision 
in the prediction of the thresholds by adding terms to the 
normal form: the minimal order of truncation of the series 
is determined by N r , that of the first resonant term in the 
normal form. However there is an optimal order that can be 
assessed by exploring the asymptotic properties of the series 
(Pucacco et al. 2008), but this issue is beyond the scope of 
the present work. In the following sections we compute the 
series (39) in the most significative cases. 
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3 BIFURCATION OF THE LOOPS 



A = M = 2, 



The system given by (1) represents motion in the symmetry 
planes of a triaxial galaxy. In each of those planes, the sym- 
metry axes directly give periodic orbits. Consider the models 
(4) at low energy: since the dynamics are slightly different 
from those of a harmonic oscillator (a = 2), we may expect 
them to be stable oscillations. What happens when energy 
increases? Nonlinear dynamics give asynchronous motions, 
with frequencies depending on amplitudes. Instability can 
be triggered by low-order resonance and we can expect a 
transition to instability and the birth of new orbits. 

The most common occurrence is that of the loops, closed 
orbits simply encircling the origin. We are going to see that 
the bifurcation providing the loops can be easily described 
with a 1:1 resonant normal form: they correspond to the 
anti-phase class tpo = ±7r introduced above (there are two of 
them, one rotating clockwise, the other counter-clockwise). 
The other class of 1:1 resonant orbits, the in-phase ?/>o = 
inclined orbits, are straight segments rotated with respect to 
the principal axes (de Zeeuw & Merritt 1983) and are forbid- 
den in the case of strictly elliptical equipotentials. Therefore 
we may ask ourselves how much we have to 'deform' the el- 
liptical equipotentials in order to accommodate for this class 
too. 



£ = Ji + J 2 , n=J 1 -.h- l V = 2(0i -0 2 ), 
the normal form (43) becomes 



(44) 



/Cn = 



By 



^TZ+ -j- (37T + (£ -7r)(2 + cos^)) , (45) 

where constant terms have been neglected for simplicity and, 
since all non-constant terms are of the same order in e, it too 
has been factored out. ICn defines a one-degree of freedom 
system with the following equations of motion 

Bi. 



n = ^-(£ 2 -^ 2 )sin^, 
4> = | + ^-71(1 - cost/;). 



(46) 
(47) 



As anticipated above, -0 = and ip = ±7r solve (46) when 
1Z 7^ ±£. However, for tf> = 0, equation (47) does not admit 
any solution in 1Z. This means that inclined orbits do not 
appear. Rather, for tp = rr, we find the solution 

n = n i = -jB-;- ( 48 ) 

In view of (44), the constraints ^ Ji , J2 ^ £ applied to this 
solution give the condition of existence for loop orbits. By 
using (40) for the ellipticity parameter, we find the threshold 



3.1 The 1:1 resonant normal form 

The general treatment of the m=n=l symmetric resonance 
with two reflection symmetries has been given by de Zeeuw 
& Merritt (1983) on the basis of previous work by Verhulst 
(1979). Their results, based on the method of averaging, in 
principle contain the answer to the questions posed at the 
start of this section. We prefer to present these results within 
the framework of Lie transform normalization because it is 
more effective in particular when studying higher-order res- 
onances. 

We approximate the frequency ratio with (16) in the 
1:1 case, 



St 



1 + 8 = 1 + 5e 2 



(40) 



so that, after the scaling transformation (11-12) and (17), 
the Hamiltonian (1) takes the form (19) with 



#o = -(pi + x{ + p 2 2 + xi) 

H 2 = 5 { x l + pl) + Bl (xl + xi)\ 



(41) 
(42) 



We truncate at order N = N r = 1 and consistently expand 
the ellipticity parameter according to (18) up to the same 
order. A standard normalization procedure (Belmonte et al. 
2007; Marchesiello & Pucacco 2011) transforms the Hamil- 
tonian (19) into the 'normal form' 



An = J1 + J2+ e 2 5J 1 + 



(43) 



( \(Jl + 4) + JiM2 + cos(20i - 2fa)) 



3.2 Bifurcation of the 1:1 resonant periodic orbits 

By introducing quasi-canonical variables adapted to the 
resonance by means of the linear combinations (34) with 



2Bi 



(49) 



To be concrete we can express this result in the case of the 
Q-models (4). In view of the rescaling and of the expansion 
of the energy as a truncated series in the parameter £, we 
have that E — lo 2 £ = S/q is a first order estimate of the 
'true' energy of the orbital motion. We can use the above 
critical values to establish the instability threshold for the 
model problem given by potentials (4): 



E^E t = 



4|T 



In the range 
0.7 < q < 1.3, 



(50) 



(51) 



which can be considered as 'realistic' for elliptical galaxies, 
the thresholds (50) give estimates correct within a 10% if 
compared to numerical computations (Belmonte et al. 2007; 
Pucacco et al. 2008). When (50) is satisfied, loop orbits bi- 
furcate from the y-axial normal mode in the oblate case, and 
from the i-axial normal mode in the prolate case (March- 
esiello & Pucacco 2011). At the same bifurcation values, the 
normal mode suffers a change of stability, passing from sta- 
ble to unstable when the new orbit is born. By direct check 
of the nature of the critical point (1Z = TZe,-^ = n) of the 
function (45) the loop, when it exist, is stable. 

To get a higher precision, we have to include higher- 
order terms in the series expansion. If we expand the poten- 
tial up to order six and truncate the normal form at iV = 2, 
the critical energy (50) up to order two in the detuning pa- 
rameter is given by 



Eu 
E21 



4 „ . 2(2 + 3a)., n2 
2 — a (a — 2) z 

4 ., , 2(5a-2)., , 2 
-2-^^-^ + -(^W^-^ 



(52) 
(53) 
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respectively in the oblate and prolate case. These generalize 
the expression for the logarithmic (a — 0) potential reported 
in (Belmonte et al. 2007). We don't give the details to ar- 
rive at these results since the second order case is explicitly 
treated in Section 4 where it is necessary to describe the 1:2 
resonance. 



3.3 Ellipse-breaking deformations 

Let us now consider a deformation of the potentials (4) by 
introducing a small parameter (3 such that 

B 1 (s 4 + 2fix 2 y 2 ), 



(54) 



with 'boxy' or 'disky' shapes of the level curves when respec- 
tively /3 < and /3 > 0. As we will show in the following, 
the presence of the parameter /3 affects the bifurcation of 
inclined orbits. The normal form of the system is the same 
as Kii in (45) except that the coefficient in front of the 
resonant term is replaced by 



Sr (1+0). 



(55) 



The important point is that the second equation of motion 
for the reduced system becomes 



(1 + j9)(2 + cobV)). 



Now, for ip = equation (56) admits the solution 
S 



(56) 



(57) 



3Bi/r 

This fixed point determines two inclined orbits for the orig- 
inal system. For -0 = n the right-hand side of equation (56) 
vanishes for 

5 



n = n t {p) = - 



Bi(2-py 



(58) 



This determines the loops as before and is only slightly 
changed with respect to the solution (48) found above. 

Working as usual in the family (4), the constraints ^ 
Ji,J 2 ^ £ translate into the existence condition 



£ > £i, 2 ;(/3) 
and 

£ 3= £iMP) 



, 4(1 -g) 



4(1 ~ g) 
(2 -«)(/? -2)' 



(59) 



(60) 



where now, with the indexes 1,2, we now distinguish be- 
tween the bifurcations from the two normal modes. The 
critical values (59) correspond respectively to the bifurca- 
tion of inclined orbits from the y, :r-axial normal mode and 
the same with (60) for the loops. This distinction is relevant 
if one is interested in which normal mode suffers a change 
of stability when a new orbit arises. 

Thus, if we break the ellipticity of the potential, in- 
clined orbits appear, however the smaller the deformation, 
the higher the threshold value (59) . Loops continue to bifur- 
cate at a lower energy: to change the bifurcation sequence, 
unreasonable high values of /? are required. The phenomenon 
is anyway interesting because it can easily be checked that 
the two families are always of different stability nature: the 
stable one is the first to appear, therefore there is a critical 
value of P at which there is an exchange of stability between 



loops and inclined. The special value j3 = 2 producing the 
singularity in (60) is associated to exact separability in ro- 
tated Cartesian coordinates which forbids the existence of 
the loops. 

One may wonder if the inclusion of additional terms in 
the series does modify qualitatively the results obtained at 
lower orders: a nice result provided by the theory of singu- 
larity (Broer et al. 1998) proves that this is not the case, at 
least for the symmetric 1:1 resonance. The case with ellipti- 
cal equipotential (/3 = 0) is in a certain sense degenerate, but 
a generic symmetry-preserving deformation is stable. The 
meaningful information is essentially contained in the nor- 
mal form truncated at N = 1 since, even adding higher-order 
terms to the original physical Hamiltonian, one can always 
find a non-linear coordinate transformation allowing us to 
eliminate the extra terms from the normal form: in other 
words the bifurcations predicted by using (45) (including 
the deformation) are qualitatively reliable and can only be 
quantitatively improved with a higher-order normalization 
(Pucacco et al. 2008). 



4 BIFURCATION OF THE BANANA AND 
ANTI-BANANA 

Another important class of bifurcations is that of banana or- 
bits (Miralda-Escude & Schwarzschild 1989) usually associ- 
ated to the instability of the major-axis orbit. It corresponds 
to a pair of in-phase (tpo — 0) oscillations with frequency ra- 
tio 1:2. The anti-phase family are the figure-eight periodic 
orbits, or anti-banana: we will show that in the potentials 
(4) stable bananas bifurcates at lower energies than unstable 
anti-bananas for relevant values of the parameters. 

In the case of the m=l,n=2 resonance with reflection 
symmetries about both axes, we know from the general ex- 
pression (32) of the normal form, that the normalization 
procedure must be pushed at least to order N r — 2. The 
terms in the series expansion (19) are now given by 



B 1 {x 2 1 +2x 2 2 f, 



H = -m(p 
H 2 = 5(x 2 4 
Hi = 25B 1 {x\-Axi)+B 2 (x 2 1 +2xlf . 
After normalization, we get the 'normal form' 



K 12 = j2 £2kK ^> 



where 



i) 



Ji +2J 2 , 



( ^Ji 2 +4J 1 J 2 +6J| , 



(61) 
(62) 
(63) 

(64) 

(65) 
(66) 



35B 1 (J 2 - 4J 2 2 ) - (17B 2 - 10B 2 ) -J't + 2J : 



1(465? - 27B 2 )JiJl - f y B 2 - 9B 2 
|(2B 2 - B 2 )JIJ 2 cos(46»! - 20 2 ). 



(67) 



We remark that in the computation of (64) and re- 
sults thereof, the use of algebraic manipulators like 
Mathematica® is practically indispensable. 
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The quasi-canonical transformation to adapted reso 
nance coordinates now is 

( Ji= £ + 211 

J 2 = 2£-TZ 

ip= 40i - 26 2 

I X = 20i + 402 

and the effective Hamiltonian 



(68) 



K. 12 {n,i>;£) = K 12 (J a {£,K),0 a (ip, X )) 
defines the following equations of motion 

n = 

ip = 2s 2 I Bi (3£ - 4TZ) - 25) + 
- ie 4 [5.4(71; £,$) - B{TZ;£) cos 



(69) 



|e 4 (273 2 - 73 2 )(2£: - TZ)(£ + 2TZ) 2 sin^, (70) 



(71) 



where 

A = 
+ 

B = 



36B 2 £(-3£ + 411) 

Bl{155£ 2 - 276£ft + 48ft 2 ) + 72Bi£S, 
9(273 2 - B 2 )(7£ 2 + 8£TZ - 12ft 2 ). 



The fixed points of this system give the periodic orbits 
of the original system. The pair of solutions with ft = 2£, 
ft — —£/2 respectively correspond to the normal modes 
along the a;-axis and j/-axis. Let us look for periodic orbits 
in general position. We start with setting ip = and looking 
for ft-solutions of ip — 0. Since we are dealing with a per- 
turbation problem in e, we look for a solution in the form 
(Henrard 1969) 



ft = fto+fti£ 2 +0(e 4 ). 



(72) 



We substitute (72) in (71) with -0 = and collect terms 
up to fourth order in e. Equating to zero the coefficient of 
second order, we find that fto has to satisfy 



73i(3£- 4ft) -25 = 
which gives 

^ 3„ 5 
1Z ^4 £ -2B[- 



(73) 



(74) 



Once computed fto we find the coefficient of the second order 
term in the expansion of the fixed point 



TZb = fto + HbiE , ip = 0, 

which determines the banana orbits: 



(75) 



Jib = £ + 2ft(,, 
J 2b — 2£ — IZb- 

Similarly, for ip — 40i 
form 



202 



TZ a = fto + TZaie , ip = ir, 



(76) 
(77) 

7r, we find a solution of the 
(78) 



and Jia = £ + 21Z a , J 2a = 2£ — TZ a , corresponding to the 
antibanana orbits. 

In view of (68), the constraints ^ Ji ^ 5£, ^ J 2 ^ 
5£ /2 applied to these solutions give the condition of exis- 
tence for these periodic orbits. Non trivial existence condi- 
tions can be found by solving Ji l2 (, ^ for the bananas and 
Ji,2a ^ for the anti-bananas. The implicit function theo- 
rem assures that there exists unique solutions £ c = £(5) in 



each cases determining the bifurcation thresholds. For the 
bananas, up to the second perturbative order we get 



_ 2 - 59gf - 275 2 ?2 2 

Sbi = -5m 5+ — T573f — S£ > 



(79) 



3673 2 J2 2 
£ , 



(80) 



2 ~ 97B 2 
£b2 = 5B; S+ 1573? 

which respectively determine the bifurcation from the x- 
axial normal mode in the first case, and from the y-axial 
normal mode in the second case (we discuss below which of 
these possibilities actually shows up). Similarly, the thresh- 
old values that gives the existence condition of anti-banana 
orbits are given by 

2 g , 19 Bf - 973 2 t2J2 
£al = -5Bl S+ 373? 5£ > (81) 

2 - 9773 2 - 3673 2 ;2 2 
£a2 = 5B7 5+ 1573? 5£ - (82) 

By comparing (80) with (82) we see a first interesting re- 
sult: if the bifurcation occur from the j/-axis, banana and 
anti-banana appear together. It is therefore important to dis- 
criminate between the two possibilities. Since the dominant 
term in the series is the first and £ must be positive, we 
see that case 1 (bifurcation from the a;-axis) or 2 (bifurca- 
tion from the j/-axis) occur if 8 and 73i have different sign 
or not. To write the expressions of the bifurcation curves 
in the physical (q, 7?)-plane, according to the rescaling (17) 
with n = 2, on the two axial orbits we have 



73i 



= 5£ £ 2 + y73i£:V + 0( e 6 ), 



E 2 = 5£e 2 + (^jBi£ 2 - 10£5j e 4 + 0(e 6 ) 

so that we get 

2 7773? - 2773 2 y2 
Ebl = -Bl 5 + W( S ' 



E a l 



2 . 11373 2 - 4573 2 

b1 5 + — w! 



for the bifurcations from the a:-axis and 
E b2 



2 J , W3Bj - 3673 2 , 2 
= E a2 = — 6+ — f S ' 



(83) 
(84) 

(85) 
(86) 

(87) 



for the bifurcations from the j/-axis. To be concrete, for our 
family (4) we have that, with a > 0, the coefficient 73i is 
negative. The ellipticity is usually q > 1/2 so that S > 0, 
therefore relevant thresholds are 



E a l = 



16 

2 - a 

16 



2 - a V 2 
Since the difference 

2 + a 



1\ 8(41q-10) / 1 
q -2j + 3(2 -a)" [ q ~2 
1\ 8(53a + 14) / 1 



3(2 -a) 2 



E a l — Eb 



32 



(2 



q -~2 



(89) 



(90) 



is positive, we verify that, for models in the class (4) and 
with parameter ranges useful for elliptical galaxies, the bifur- 
cation sequence is always from the major axis, with bananas 
appearing at lower energies than anti-bananas. By checking 
the nature of the two critical points (75,78), it can be seen (it 
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is a tedious but straightforward computation, Marchesiello 
& Pucacco (2012)) that in systems (4) the first family is al- 
ways stable and the second unstable: (88,89) generalize the 
corresponding expressions for the logarithmic (a = 0) po- 
tential reported in Belmonte et al. (2007) . As long as the ba- 
nana does not bifurcate the major axis is stable and parents 
'box' orbits. It loses its stability at the first bifurcation and 
regains it at the second. It is natural to ask how much these 
results are affected by ellipse-breaking deformations: we can 
say that, in analogy with what seen for the 1:1 resonance, 
the hierarchy of bifurcations changes only for unreasonable 
high values of the deformation parameter. 



lution of the form 



5 HIGHER-ORDER SYMMETRIC 
RESONANCES 

As is well known (Miralda-Escude & Schwarzschild 1989) 
stable periodic orbits corresponding to higher-order reso- 
nances and quasi-periodic orbits parented by them give a 
small but not-negligible contribution to regular dynamics 
in systems with cores. In realistic cases with mixed (regu- 
lar+chaotic) dynamics it is conjectured that these 'boxlets' 
may become important in shaping the bulk of the density 
distribution (Zhao 1999; Zhao et al. 1999). The main differ- 
ence of these families from those seen above consists of the 
fact that their bifurcation is not connected with the loss (or 
regain in case of a second bifurcation) of stability of the nor- 
mal mode. The birth of periodic orbits with N r > 2 is rather 
due to breaking of a resonant torus around the normal mode 
and is correctly described by applying the Poincare-Birkhoff 
theorem (Arnold 1989): however, the technique we applied 
above continues to work and the conditions for the exis- 
tence and stability of an m/n-resonant periodic orbit with 
m + n > 3 can still be found by constructing the appro- 
priate normal form and locating fixed points of the reduced 
system. 

A technical issue worth to be clarified is the following: 
by reducing the resonant normal form (32) truncated at or- 
der N r by means of the transformation (34), we obtain a 
polynomial of degree N r + 1 in 1Z. The corresponding equa- 
tion of motion for i/j produces a pair of algebraic equations 
of degree N r which have to be solved to locate the fixed 
points (one for each solutions tpo, cfr. point ii in Sect. 2). 
This problem is very difficult to solve if, for N r > 2, we 
aim at general solutions depending on the parameters of the 
system. However we are not interested in every solution but 
only in those connected with the passage through the cho- 
sen resonance. We can therefore resort to the perturbation 
method we have described in detail in the previous section 
on the 1:2 resonance. In that case, with iV r = 2 we had to 
solve two equations of second degree (cfr. the rhs side of 
(71)): this clearly does not represent a problem since we can 
write explicitly the two pair of solutions. However, in each 
pair, only one solution is geometrically acceptable because 
it satisfies the condition at resonance; the other must be 
discarded by direct check. The perturbative method based 
on the construction of the series (72) (Henrard 1969) au- 
tomatically selects the acceptable solution. The method is 
therefore extremely useful for higher-order resonances: a so- 



N r -1 



TZ k e +U(e 



(91) 



easily allows us to select the meaningful solution without 
any loss in accuracy. 

We have applied the method to the case of fish orbits 
corresponding to the (anti-phase) 2:3 resonance. In this case, 
N r — 4: the Hamiltonian series must be expanded up to in- 
clude terms of degree 10 (B4 in the original potential). The 
explicit expressions of the normal form in the general class 
(5) and for the family (4) are a bit heavy to write and are 
reported elsewhere (Marchesiello 2012): they are available 
upon request as Mathematica® notebooks. Anyway the pro- 
cedure is a straightforward extension of that illustrated in 
the previous section. 

The threshold for the existence of fish orbits turns out 
to be 



E f = 



3 9 
2Bi 80S 3 



(149S? - 60B 2 )S 2 



27(76715* - 7840B 2 B 2 + 3600B| - 1500BiB 3 )<5 3 



81 



1600Bf 

r (4852431Bi - 8889450BiB 2 



448000B! 7 

+9116400B 2 B 2 - 3780000S| - 36260005? B 3 
+315OOOOB1B2B3 - 4900005? B 4 )<5 4 . (92) 



This result is undoubtedly unpleasant to write (and read!) 
but it testifies what is the rule with high-order expansions. 
However, trusting the normalization program and paying 
attention to write down the results without errors, the series 
give us numbers we can use in specific cases. In terms of the 
parameters of the family (4) we get 



Ej 



12 



9(22 + 69a) 2 



+ 



2 - a 10(2 - 
9(4372 + 2508a + 4853a 2 )S 3 
200(2 - a) 3 

27(1368856 + 3109116a + 542642a 2 + 1468293a 3 )5 4 



56000(2 



(93) 



where in this case 



6 = q 



(94) 



This result completes and generalizes the treatment of the 
logarithmic case presented in Belmonte et al. (2007). We 
may ask if it is worth the effort: in the logarithmic case 
(a = 0), Miralda-Escude & Schwarzschild (1989) numeri- 
cally found E f (q = 0.7) = 0.21 and E f {q = 0.9) = 2.28 
that we can consider experimental exact threshold values. 
Our analytic result predicts Ef(q = 0.7) = 0.206 and 
Ef(q — 0.9) = 2.10. The agreement is excellent near the 
resonance (S — 0.7 — 2/3 ~ 0.03) and only moderate further 
away from it (S = 0.9 — 2/3 ~ 0.23). However, we remark 
that the energy level E = 2.28 is extremely high if seen with 
the eye of the perturbation theorist: an error of 8% may then 
appear not so bad. Moreover it is possible to improve the 
quality of the prediction by going to still higher orders. 

If one is only interested in a rough prediction around 
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a general m:n resonance (Pucacco 2009), from these results 
we can deduce the general first order expression 

E m:n = -^—8, (95) 
mBi 

that, for the family (4), gives 

E m -. n = *» (96) 
m(2 — a) V n ) 

The example of the 3:4 resonance (the pretzel) is a 
good test: for the logarithmic potential Miralda-Escude & 
Schwarzschild (1989) numerically found E Z: 4,(q — 0.7) = 
0.25 and E z . A {q = 0.9) = 1.22. Eq.(96) with a = pre- 
dicts £ 3 :4(q = 0.7) = 0.27 and E a - A (q = 0.9) = 0.80. The 
agreement is quite good near the resonance; moreover, since 
for q < 3/4 (96) is negative, accordingly with the treat- 
ment of the previous cases, we may predict that in the case 
q — 0.7 the 'bifurcation' is from the t/-axis, as actually found 
by Miralda-Escude & Schwarzschild (1989). 



6 DISCUSSION 

In the present section we discuss some implications of the 
results described above and present open problems and pos- 
sible directions to cope with them. Here we also recall that 
the approach we have followed is not the only possible and 
that both the normalization and the reduction can be ob- 
tained by exploiting alternatives like the Lissajous trans- 
formation (Deprit & Elipe 1991), the method of geometric 
invariants (Hanfimann & Sommer 2001) and the singularity 
theory (Broer et al. 1998) mentioned in the introduction. 

6.1 Asymptotic expansions 

Series like those described in this work are asymptotic: this 
means that a truncation of the series, say at order N, appar- 
ently converges in a given domain only for iV < N op t, the 
optimal truncation order linked to the extent of the domain. 
We remark that this semi-convergence is in general not as- 
sociated to a true function: rather, it is only associated to 
a local geometric object we use as an invariant surface in 
the regular part of phase space. The optimal truncation de- 
pends on the problem at hand and to assess it a priori is 
quite difficult (Efthymiopoulos et al. 2004). In Pucacco et 
al. (2008) we have tried to estimate A?opt for two members 
of the family (4), a — 0,1. For the bifurcation of the ba- 
nana in the logarithmic potential, we obtained N op t > 7 for 
q < 0.7, iVopt = 6 for q = 0.8 and jV opt = 3 for q = 0.9. 
In this case (the worst being the furthest from exact reso- 
nance) the relative error of the prediction is 11%. However, 
the quality of the prediction (and the corresponding opti- 
mal order) can be further improved if different techniques of 
summation are employed. Scuflaire (1995) suggested to use 
the continued fraction method (Bender & Orszag 1978) to 
re-sum asymptotic series: we applied this idea to the bifur- 
cation threshold series and in the worst case just mentioned 
(banana with a — 0, q = 0.9) we got N op t = 5 lowering the 
relative error to less than 4%. What is indeed remarkable in 
this result is that the bifurcation energy is E = 3.6. For the 
logarithmic potential this corresponds to a radius of order 
40 times larger than the convergence radius of the original 



series (5) so that we have an outstanding evidence of the 
power of asymptotic expansions. 

6.2 Stackel fits 

Separable systems play an important role among integrable 
systems since they provide explicit solutions for the orbit 
structure. The application of Stackel systems to approxi- 
mate the dynamics of galaxies is therefore a classical field 
(van de Hulst 1962; de Zeeuw 1985b; Kent & de Zeeuw 1991; 
van de Ven et al. 2008). de Zeeuw & Lynden-Bell (1985) pro- 
posed a 'Stackel fit' of galactic potentials around an equi- 
librium to take advantage of the opportunity of exploiting 
the integrals of motion of systems separable in elliptical co- 
ordinates. At order N = 1 the number of free parameters is 
sufficient to fit any expansion; at higher orders the fit is con- 
strained by conditions on the coefficients. The method works 
since the dynamics of a Stackel system separable in elliptical 
coordinates resemble that of the 1:1 resonance for potentials 
of the form (5): however, the results obtained in Subsection 

3.3 warn us from excessive confidence in the method. In fact 
we can fit a potential of the form (54) or even more general 
in which inclined orbits may play a role: however the fitting 
Stackel potential does not support inclined. Although this 
may not be of particular relevance in galactic applications, 
it is a problem as a matter of principle. We remark that 
Stackel systems do not end with those mentioned above, but 
include those separable in other coordinate systems. In 2 di- 
mensions, separability in parabolic coordinates can be used 
to model elliptical disks (Sridhar & Touma 1996, 1999): in 
this case there is a relation with the 1:2 resonance. How- 
ever, systems separable in parabolic coordinates accommo- 
date bananas and quasi periodic orbits parented by them, 
but do not support their anti-phase companions. 

6.3 Surfaces of section 

By inverting the transformation leading to the normal form 
we can compute formal integrals of motion (Contopoulos et 
al. 2003; Contopoulos 2004) which have to be interpreted as 
asymptotic series as prescribed in Subsection 6.1. The most 
immediate use of these expansions is to construct approx- 
imations of Poincare surfaces of sections: for the logarith- 
mic potential, Belmonte et al. (2007) show that, at suffi- 
ciently high energy, surfaces constructed around low-order 
resonances display a quite close resemblance with those nu- 
merically obtained in the scale-free limit by Miralda-Escude 
& Schwarzschild (1989). Moreover, by using asymptotic se- 
ries as true phase-space conserved functions in a suitable 
domain, bifurcation curves can be computed by investigat- 
ing the nature of the critical points of these functions. The 
results are identical to those obtained with the normal form 
when expressed as series in the detuning: either approaches 
being effective, one can chose which minimize the computa- 
tional effort. 

6.4 Order and chaos 

The domain of 'semi-convergence' of asymptotic series ap- 
proximating invariant surfaces of generic systems can be 
taken as a measure of their regular dynamics. We have seen 
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that, as a matter of principle, regular phase-space zones 
associated to resonances of any order can be adequately 
included and described. The approach to high-order reso- 
nances is dual: either their role is considered to be marginal 
(Sanders 1976) or they are considered as an inescapable 
signature of chaos (Binney & Tremaine 2008). However, in 
several interesting cases (see e.g. the scale-free models with 
a > treated by Touma & Tremaine (1997)) we have that 
different resonances coexist without overlapping for a large 
range of parameters. Resonance manifolds generate a struc- 
ture that can be understood via reduction (Tuwankotta & 
Verhulst 2000). Regular dynamics are 'complicated' but def- 
initely not chaotic, so efficient tools to investigate their fea- 
tures are extremely useful. 



6.5 3D models 

The most relevant generalization is towards 3 dimensional 
systems. The pioneering work by de Zeeuw (1985a) still re- 
mains a major contribution since mathematicians, although 
have devoted much effort to this issue, analyzed in general 
only simple abstract models (Sanders et al. 2007). de Zeeuw 
(1985a) gave an almost complete study of the orbit structure 
of a generic quartic potential around the 1:1:1 resonance. 
The relevance of this case is testified by the fact that, in 
spite of a radical change in our understanding of elliptical 
galaxies with cusps affecting their overall dynamics, the two 
orbit families characterizing triaxial systems are still con- 
sidered to be the boxes and the long axis-tubes (van den 
Bosch & de Zeeuw 2010): we therefore see that the study 
of the stability of the s-normal mode and the condition for 
existence of stable loops in the yz-plane as studied in this 
work is very useful. 

The main problem with 3 degrees of freedom is that the 
normal form itself is in general not integrable: the normaliza- 
tion procedure of resonant Hamiltonians provides only one 
formal integral (Gustavson 1966) in addition to energy. How- 
ever, the study of the stability of the three normal modes 
and the bifurcations of periodic orbits in general position 
can be done even in the absence of a third integral. The 
step towards a general analysis of relevant cases like the 
1:2:2 and 1:2:3 resonances seems to be within the reach of 
the method. We also recall that a small bulk rotation of the 
ellipsoid can be included with a suitable canonical transfor- 
mation (de Zeeuw & Merritt 1983). 



7 CONCLUSIONS 

We have presented a general analysis of the orbit structure 
of 2D potentials with self-similar elliptical equipotentials. 
The main results are the following: 

The 1:1 resonance is associated to the appearance of the 
loops and leads to the destabilization of the y-axis orbit in 
the oblate case and of the rr-axis orbit in the prolate case. 
Inclined orbits are never present and may appear only when 
the equipotentials are heavily deformed. 

The 1:2 resonance determines the appearance of ba- 
nanas and anti-banana orbits: the first family is stable and 
always appears at a lower energy than the second, which is 
unstable. The bifurcation sequence produces the change in 



the stability character of the major axis orbit and is modified 
only by very large deformations of the equipotentials. 

Higher-order resonances appear at intermediate ener- 
gies which can be predicted with good accuracy. 

We have analyzed several issues connected with the ap- 
proach and sketched the directions for further work. In par- 
ticular, we think that evaluating the overall predictive power 
of the method based on asymptotic expansions is a decisive 
step if one is interested in studying stationary or rotating 
triaxial potentials. 
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